%     VIRTUAL FABRIC - программа для получения эффективных свойств текстильных композитов
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function Write_Input_Thermal_Yarn(Job_Name, Work_Directory, NODES, ELEMENTS, Material_Properties, Coordinate_Systems, RVE_Coordinates, Node_Merge_Tolerance)

global ErrCode;

node_coord_format = '%.12e';
orientation_format = '%.12e';
mat_prop_format = '%.12e';
equation_format = '%.12e';
bc_format = '%.12e';

Part_Name = 'UD_RVE';
Model_Name = Job_Name;

% Проверка нумерации узлов
for i=1:size(NODES,1)
    if (NODES(i,1) ~= i)
        fprintf (1, '\nОшибка функции Write_Input_Thermal_Yarn: некорректная нумерация узлов!\n');
        ErrCode = -1;
        return;
    end
end

% Проверка нумерации элементов
for i=1:size(ELEMENTS,1)
    if (ELEMENTS(i,1) ~= i)
        fprintf (1, '\nОшибка функции Write_Input_Thermal_Yarn: некорректная нумерация элементов!\n');
        ErrCode = -1;
        return;
    end
end

File_name = Job_Name;
File_path = [Work_Directory '\' File_name '.inp'];

Fout = fopen(File_path, 'wt');

if (Fout < 0)
    fprintf (1, '\nОшибка функции Write_Input_Thermal_Yarn: не удалось открыть для записи файл [%s] !\n', File_path);
    ErrCode = -1;
    return;
end

fprintf (Fout, '*Heading\n');

fprintf (Fout, '** Job name: %s Model name: %s\n', Job_Name, Model_Name);
fprintf (Fout, '** Generated by: Virtual_Fabric\n');

fprintf (Fout, '*Preprint, echo=NO, model=NO, history=NO, contact=NO\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** PARTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Part, name=%s\n', Part_Name);

fprintf (Fout, '*Node\n');

for i=1:size(NODES, 1)
    fprintf (Fout, ['%i, ' node_coord_format ', ' node_coord_format ', ' node_coord_format '\n'], NODES(i,1), NODES(i,2), NODES(i,3), NODES(i,4));
end

DC3D8_Set = zeros(size(ELEMENTS,1),1);
DC3D20_Set = zeros(size(ELEMENTS,1),1);

counters = zeros(2,1);

for i=1:size(ELEMENTS,1)
    
    switch ELEMENTS(i,3)

        case 8

            counters(1) = counters(1) + 1;
            DC3D8_Set(counters(1)) = i;

        case 20

            counters(2) = counters(2) + 1;
            DC3D20_Set(counters(2)) = i;

    end
    
end

if (counters(1) > 0)
    DC3D8_Set = DC3D8_Set(1:counters(1));
else
    DC3D8_Set = [];
end

if (counters(2) > 0)
    DC3D20_Set = DC3D20_Set(1:counters(2));
else
    DC3D20_Set = [];
end

if (~isempty(DC3D8_Set))
    
    fprintf (Fout, '*Element, type=DC3D8\n');
    
    for i=1:size(DC3D8_Set,1)
        fprintf (Fout, '%i, %i, %i, %i, %i, %i, %i, %i, %i\n', ELEMENTS(DC3D8_Set(i),1), ELEMENTS(DC3D8_Set(i),5), ELEMENTS(DC3D8_Set(i),6), ELEMENTS(DC3D8_Set(i),7), ...
                                                               ELEMENTS(DC3D8_Set(i),8), ELEMENTS(DC3D8_Set(i),9), ELEMENTS(DC3D8_Set(i),10), ELEMENTS(DC3D8_Set(i),11), ...
                                                               ELEMENTS(DC3D8_Set(i),12));
    end
    
end

if (~isempty(DC3D20_Set))
    
    fprintf (Fout, '*Element, type=DC3D20\n');
    
    for i=1:size(DC3D20_Set,1)
        fprintf (Fout, '%i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i,\n', ELEMENTS(DC3D20_Set(i),1), ELEMENTS(DC3D20_Set(i),5), ELEMENTS(DC3D20_Set(i),6), ELEMENTS(DC3D20_Set(i),7), ...
                                                                                           ELEMENTS(DC3D20_Set(i),8), ELEMENTS(DC3D20_Set(i),9), ELEMENTS(DC3D20_Set(i),10), ELEMENTS(DC3D20_Set(i),11), ...
                                                                                           ELEMENTS(DC3D20_Set(i),12), ELEMENTS(DC3D20_Set(i),13), ELEMENTS(DC3D20_Set(i),14), ELEMENTS(DC3D20_Set(i),15), ...
                                                                                           ELEMENTS(DC3D20_Set(i),16), ELEMENTS(DC3D20_Set(i),17), ELEMENTS(DC3D20_Set(i),18), ELEMENTS(DC3D20_Set(i),19));
        fprintf (Fout, '   %i, %i, %i, %i, %i\n', ELEMENTS(DC3D20_Set(i),20), ELEMENTS(DC3D20_Set(i),21), ELEMENTS(DC3D20_Set(i),22), ELEMENTS(DC3D20_Set(i),23), ELEMENTS(DC3D20_Set(i),24));
    end
    
end

Fiber_Elements = zeros(size(ELEMENTS,1),1);
Matrix_Elements = zeros(size(ELEMENTS,1),1);

counters = zeros(2,1);

for i=1:size(ELEMENTS,1)
    
    switch ELEMENTS(i,4)
        
        case 1
            
            counters(1) = counters(1) + 1;
            Fiber_Elements(counters(1)) = i;
            
        case 2
            
            counters(2) = counters(2) + 1;
            Matrix_Elements(counters(2)) = i;
        
    end
    
end

if (counters(1) > 0)
    Fiber_Elements = Fiber_Elements(1:counters(1));
else
    Fiber_Elements = [];
end

if (counters(2) > 0)
    Matrix_Elements = Matrix_Elements(1:counters(2));
else
    Matrix_Elements = [];
end

fprintf (Fout, '*Nset, nset=ALL_NODES, generate\n');
fprintf (Fout, ' %i, %i, 1\n', 1, NODES(size(NODES,1),1));

fprintf (Fout, '*Elset, elset=ALL_ELEMENTS, generate\n');
fprintf (Fout, ' %i, %i, 1\n', 1, ELEMENTS(size(ELEMENTS,1),1));

if (~isempty(Fiber_Elements))
    fprintf (Fout, '*Elset, elset=Fiber_Elements\n');
    Write_Set(Fout, Fiber_Elements);
end

if (~isempty(Matrix_Elements))
    fprintf (Fout, '*Elset, elset=Matrix_Elements\n');
    Write_Set(Fout, Matrix_Elements);
end

format_line = [orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format '\n'];

fprintf (Fout, '*Orientation, name=Ori_Coord\n');
fprintf (Fout, format_line, Coordinate_Systems{1}.Point_On_Axis_1(1), Coordinate_Systems{1}.Point_On_Axis_1(2), Coordinate_Systems{1}.Point_On_Axis_1(3), ...
                            Coordinate_Systems{1}.Point_On_Plane_12(1), Coordinate_Systems{1}.Point_On_Plane_12(2), Coordinate_Systems{1}.Point_On_Plane_12(3), ...
                            Coordinate_Systems{1}.Origin(1), Coordinate_Systems{1}.Origin(2), Coordinate_Systems{1}.Origin(3));
fprintf (Fout, '3, 0.\n');

fprintf (Fout, '** Section: Fiber_Sectior\n');
fprintf (Fout, '*Solid Section, elset=Fiber_Elements, orientation=Ori_Coord, material=Fiber_Material\n');
fprintf (Fout, ',\n');

fprintf (Fout, '** Section: Matrix_Section\n');
fprintf (Fout, '*Solid Section, elset=Matrix_Elements, material=Matrix_Material\n');
fprintf (Fout, ',\n');

fprintf (Fout, '*End Part\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** ASSEMBLY\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Assembly, name=Assembly\n');
fprintf (Fout, '**\n');
fprintf (Fout, '*Instance, name=%s-1, part=%s\n', Part_Name, Part_Name);
fprintf (Fout, '*End Instance\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node\n');
fprintf (Fout, '      1,           0.,           0.,           0.\n');
fprintf (Fout, '      2,           0.,           0.,           0.\n');
fprintf (Fout, '      3,           0.,           0.,           0.\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_x\n');
fprintf (Fout, ' 1,\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_y\n');
fprintf (Fout, ' 2,\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_z\n');
fprintf (Fout, ' 3,\n');

RVE_AG_NegX_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,1)+Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);
RVE_AG_PosX_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(1,2)+Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

RVE_AG_NegY_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,1)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);
RVE_AG_PosY_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

RVE_AG_NegZ_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,1)+Node_Merge_Tolerance);
RVE_AG_PosZ_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(3,2)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

if (isempty(RVE_AG_NegX_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegX_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_NegY_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegY_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_NegZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegZ_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosX_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosX_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosY_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosY_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosZ_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegX_RVE_AGSORT) ~= numel(RVE_AG_PosX_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegX_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosX_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegY_RVE_AGSORT) ~= numel(RVE_AG_PosY_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegY_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosY_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegZ_RVE_AGSORT) ~= numel(RVE_AG_PosZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegZ_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosZ_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

RVE_AG_PosX_RVE_AGSORT = Match_Nodes(RVE_AG_PosX_RVE_AGSORT, RVE_AG_NegX_RVE_AGSORT, NODES, 1, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

RVE_AG_PosY_RVE_AGSORT = Match_Nodes(RVE_AG_PosY_RVE_AGSORT, RVE_AG_NegY_RVE_AGSORT, NODES, 2, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

RVE_AG_PosZ_RVE_AGSORT = Match_Nodes(RVE_AG_PosZ_RVE_AGSORT, RVE_AG_NegZ_RVE_AGSORT, NODES, 3, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

fprintf (Fout, '*NSET, NSET=RVE_AG_NegX_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegX_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosX_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosX_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_NegY_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegY_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosY_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosY_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_NegZ_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegZ_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosZ_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosZ_RVE_AGSORT);

fprintf (Fout, '**\n');
fprintf (Fout, '** EQUATIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Equation\n');

fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosX_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegX_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_x, 11, ' equation_format '\n'], RVE_Coordinates(1,2) - RVE_Coordinates(1,1));
fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosY_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegY_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_y, 11, ' equation_format '\n'], RVE_Coordinates(2,2) - RVE_Coordinates(2,1));
fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosZ_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegZ_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_z, 11, ' equation_format '\n'], RVE_Coordinates(3,2) - RVE_Coordinates(3,1));

fprintf (Fout, '*End Assembly\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** MATERIALS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Material, name=Fiber_Material\n');
fprintf (Fout, '*Conductivity, type=ORTHO\n');
fprintf (Fout, [mat_prop_format ', ' mat_prop_format ', ' mat_prop_format '\n'], Material_Properties.Fiber_lambda_11, Material_Properties.Fiber_lambda_22, Material_Properties.Fiber_lambda_33);
fprintf (Fout, '*Density\n');
fprintf (Fout, [mat_prop_format ',\n'], Material_Properties.Fiber_rho);
fprintf (Fout, '*Specific Heat\n');
fprintf (Fout, [mat_prop_format ',\n'], Material_Properties.Fiber_Cp);

fprintf (Fout, '*Material, name=Matrix_Material\n');
fprintf (Fout, '*Conductivity\n');
fprintf (Fout, [mat_prop_format ',\n'], Material_Properties.Matrix_lambda);
fprintf (Fout, '*Density\n');
fprintf (Fout, [mat_prop_format ',\n'], Material_Properties.Matrix_rho);
fprintf (Fout, '*Specific Heat\n');
fprintf (Fout, [mat_prop_format ',\n'], Material_Properties.Matrix_Cp);

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_X-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_X-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, ['RP-ThermGrad_x, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(1,2) - RVE_Coordinates(1,1)));

fprintf (Fout, '** Name: HomogK_X-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, 'RP-ThermGrad_y, 11, 11\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, 'RP-ThermGrad_z, 11, 11\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_Y-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_Y-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_x, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, ['RP-ThermGrad_y, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(2,2) - RVE_Coordinates(2,1)));

fprintf (Fout, '** Name: HomogK_Y-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_z, 11, 11\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_Z-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_Z-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_x, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_y, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, ['RP-ThermGrad_z, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(3,2) - RVE_Coordinates(3,1)));

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fclose(Fout);
